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ABSTRACT 

We introduce a new statistic u>i(r s ) for measuring and analyzing large-scale structure and partic- 
ularly the baryon acoustic oscillations. u>e(r s ) is a band-filtered, configuration space statistic that is 
easily implemented and has advantages over the traditional power spectrum and correlation function 
estimators. Unlike these estimators, U£(r s ) can localize most of the acoustic information into a single 
dip at the acoustic scale while also avoiding sensitivity to the poorly constrained large scale power 
(i.e., the integral constraint) through the use of a localized and compensated filter. It is also sensitive 
to anisotropic clustering through pair counting and does not require any binning of data. We mea- 
sure the shift in the acoustic peak due to nonlinear effects using the monopole u>o(r s ) derived from 
subsampled dark matter catalogues as well as from mock galaxy catalogues created via halo occupa- 
tion distribution (HOD) modeling. All of these are drawn from 44 realizations of 1024 3 particle dark 
matter simulations in a l/i -1 Gpc box at z=l. We compare these shifts with those obtained from the 
power spectrum and conclude that the results agree. We therefore expect that distance measurements 
obtained from oJo(r s ) and P(k) will be consistent with each other. We also show that it is possible to 
extract the same amount of acoustic information by fitting over a finite range using either wo(r s ) or 
P(k) derived from equal volume surveys. 

Subject headings: distance scale — cosmological parameters — large-scale structure of universe — 
cosmology: theory 



1. INTRODUCTION 

Baryon acoustic oscillations (BAOs) are relic imprints 
on the baryon distribution left behind by the interaction 
between radiation and matter in the primordial universe. 
The large amount of radiation pressure due to high tem- 
peratures in the early universe pushed matter apart as 
it congregated under the influence of gravity. Oscillat- 
ing sound waves were set up in th e primordial plasma 



due to these two com peting e ffects (IPeebles fc Yul 
Sunvaev fe Zeldovichl [T970t iBond fe Efstathiou 



1970; 



1984; 



Holtzmanl Il989t IHu fc Sugivamal 1199ft IHu fc Whitel 
19961; lEisenstein fc Hul 119981) . The distance traveled 
by these sound waves before recombination is about 
150 comoving Mpc and is known as the acoustic 
scale. When the photons stream off, the baryons 
are deposited at these characteristic separations and 
hence we would expect excesses with this separa- 
tion in the baryon distribution today. This charac- 
teristic scale makes for a very useful standard ruler 
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Blake & Glazebrookl 12001 IHu & Haiman 


2003; 


Linder 


12003: Seo & Eisenstein 2003: Matsubaral 


2004; 


Amendola et al.ll2005lL 



As the universe evolves, primordial density per- 
turbations begin to grow nonlinearly, especially on 
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small s cales (|Jain fc Bertschingeil 119941: iMeiksin et al 



1999; 



1999) 



IMeiksin fc Whitel 
This causes the 



19991 ; iScoccimarro et al 



causes tne acoustic scale to become 
slightly shifted from its pred ic ted linear theory po- 



sition (ISeo fc Eisenstein! 120051: IHuff et al.l l2007t 



20071: lAneulo et al l 120081: ICrocce fc Scoccimarro 
Sanchez et al.l 120081: ISeo et all 120081: iSmith et al 



200c 



200c 



Padmanabhan & White 200^ One must calibrate this 



shift before the acoustic scale can be used as a high 
precision standard ruler. 

Previous calibrations and other analyses of the acous- 
tic signature have been conducted using the traditional 
power spectrum and correlation function methods. Al- 
though analytically they are both perfectly adequate, the 
estimators used to derive them from simulations and ob- 
servational data are subject to numerous difficulties. 

The largest survey scales are always poorly con- 
strained due to effects such as the integral constraint 
making it a challenge to estimate the correlation func- 
tion £(r) accurately at these scales. The integral con- 
straint arises due to the fact that we do not know 
the cosmic number density of any population of mass 



tracers (Ide Lapparent et al. 


H988HBaumeart & Frvlll991l: 


IPeacock & Nicholson! 11991 


; Hamilton! 1993). Many 



1973 



1988 



techniques used to estimate £(r) (ireebiesi 
Sharol 11971: IHewettl fl98l IBlanchard fc Aliminil „ 
Landv fc Szalavl 11991 lHamiltonlll993lf take the number 
density of tracers in the survey volume to be the true 
number density. This assumption effectively ignores all 
power at scales larger than the survey size while simulta- 
neously increasing the correlation between scales smaller 
than the survey size which causes the off-diagonal co- 
variance matrix terms to be larger than they would be 
otherwise. 



2 



Limited survey volume and awkward survey bound- 
aries are the major concerns when trying to es- 
timate the power spectrum P(k). Typically, the 
measured power spectrum is a convolution of the 
window function, the Fourier transform of the se- 



true 1 


jower spectrum ( 


Feldman et all 11994: iPark et al.l 


19921 


Baumeart & Frvl 


199lHPeacock & Nicholsonlll991t 


Kaiser & Peacocklll991 


). Therefore, these P(k) estima- 



tors are biased. In the limit of infinite volume, the win- 
dow function should be a delta function. However, real 
surveys have finite volume and hence the window func- 
tion has a finite albeit very small width. This induces an 
artificial smoothing at small separations in k when at- 
tempting to deconvolve the window function from the ob- 
served density field. Pair counting estimators, like those 
for £(r), avoid this issue because the relative positions 
of all tracer particles and hence the true distribution is 
recorded. 

Estimating £(r) and P(k) also requires the binning of 
data. When any binning process is invoked, one must 
carefully define any averaging used to derive the value 
in each bin. To reduce these errors, bin sizes are made 
smaller which increases the dimensionality of the covari- 
ance matrix, making it even more difficult to estimate. 
Also, survey boundaries need to be addressed with spe- 
cial care. 

In the study of large scale structure, we are also in- 
terested in any anisotropics in the distribution of ob- 
j ects. These can result from the bulk motions of obj ects 
(IPeeblesI [l980t IDavis fe Peeblesl [1981 IKaiserl Il987h as 
well as assumed models for the Hubble parameter H{z) 
and the angular diameter distance Da(z) while calculat- 
ing object separations along the line of sight (LOS) and 
transverse directions respectively. Since, at large scales 
we would expect the distribution to be isotropic, any 
anisotropy can be deconstructed into velocity field infor- 
mation which in turn provides us with constraints on f2 m , 
the derivative of the growth function / = dlnD / d\na and 
the a nisotropic parameter = f / b, where b is the g alaxy 



tne anisotropic paramet er P = J lb, wnere b is tnc galaxy 
bias (IKaiserl 119871: iHamiltonl 119981: iPeacock et al.M200li: 
Hawkins et al.ll2003t iTegmark et al.l 120021 iZehayi et al.1 



2002t ISong fe Percivall 120081: IPercival fe White! 120091 
White et al.l 12009(1 After modeling these anisotropies, 



any residua l anisotropies can be us ed to constrain 
D A {z)H(z) (jAlcock fc Paczvnskil[l979T) . If our assumed 
H(z) or Da(z) models are incorrect then the distribution 
of objects will still appear anisotropic after the bulk mo- 
tion effects are removed. The magnitude of this residual 
anisotro py can be used to infer the true underlying cos- 
molo gy (|Padmanabhan fc White! I2008t lOkumura et al.1 
l2008h . 

Anisotropic information can be extracted from the 
redshift-space correlation function and to a lesser extent 
from the redshift space power spectrum ([Szalav et al.l 
1998). If one imagines a wide angle survey, P{k) estima- 
tors that rely on a Fourier transform from a Cartesian 
grid will suffer because an arbitrary wa ve vector (fc) will 
not necessarily be parallel to the LOS (ICole et al.l[l994 
[19951: IZaroubi fc Hoffmanl [19951 ISzaoudil 12004^ This 

means that each fc mode will contain information about 
both the LOS and transverse distributions. After aver- 
aging to obtain the spherically averaged power spectrum, 



any anistropies in the distribution of survey objects will 
have been erased. Instead of taking a spherical average, 
one can also analyze the full 3D P(k) through Legen- 
dre decomposition into radial and angular components. 
However, an infinite sum is re quired and applicability 
is limited to the lin ear regime ([Heavens fc Tavlorl [19951 : 
IPercival et al.l 120041) . Pair count estimators for £(r), on 
the other hand, record the distribution of tracers ac- 
curately because they record each pair's angle to the 
LOS as well as their separation. This means that any 
anisotropies in the distribution of galaxies will become 
obvious. Statistics sensitive to anisotropic clustering arc 
desirable as they offer us a means to probe the underlying 
cosmology. 

It is also aesthetically pleasing, to localize the acoustic 
information into a single feature at the acoustic scale. 
This is true for £(r) but not for P(k), which has oscilla- 
tory acoustic features. 

Although all of the above mentioned disadvantages of 
£(r) and P(k) are minor, it is still beneficial to derive a 
new statistic that does away with as many of the above 
setbacks as possible. w^(r s ) is an example of such an 
alternative. 

The organization of this paper will be as follows. In Sj2j 
we introduce ujg(r s ), including its properties and compu- 
tation. In [21 we describe the simulations, halo occupa- 
tion models and analysis methods we use to derive acous- 
tic peak shifts through implementation of the monopolc 
wo (r s ) . We demonstrate the mutual consistency between 
the peak shifts measured from the same simulations using 
uJo(r s ) and P(k) in ^3) This is indicative of the agree- 
ment we expect between distance measures from uo(r s ) 
and P(fc). We also show that with a reasonable finite fit- 
ting range and our current choice of filter for computing 
Wf(r s ), we can extract the same amount of acoustic in- 
formation using either ujo(r s ) or P{k) from equal volume 
surveys. Finally, we state the main results of this paper 
in Sj5] Details of the matter and galaxy model P{k) re- 
sults discus sed in this paper can b e found in the compan - 
lon papers iSeo et all ([20101 ) and iMehta et all (|in prepl ) 
respectively (hereafter referred to as S&M collectively). 

2. THE 0Ji(r s ) STATISTIC 

2.1. Equations and Properties 

We expand the angle dependence of the power spec- 
trum and correlation function out as a series of Lcgendrc 
polynomials in fj, = f ■ z = cos(#), where 9 is the LOS 
angle: 



£(r,/i) = £>(r)L/(A0 



(1) 



A 2 (fc, M )^^i^=^A^(fc)L^) (2) 



2tt 2 



so that 



f A?(fc)j,(fcr) 



(3) 



where ji is the spherical Bessel function of order £ and 
L# is the Legendre polynomial of order I. 

Imagine we have a filter, Wg(r, [i, r s ) = Wg(r, r s )Li(^,), 
which wc take to be compact and compensated 
(/ r 2 dr We(r, r s ) — 0) with a characteristic scale r s . Wc 



New Statistic for Analyzing BAO 



3 



define our statistic as the redshift-space correlation func- 
tion, £ s (r, fi) , convolved with the filter as a function of 
filtering scale r s . 



UJi{r s ) = i e / d 3 r £ s (r, n)Wi(r,r s )Lt(n) 



r 2 dr &(r)We(r, r s ) 



2£ + l 



dk 



Aj{k)W e (k,', 



(4) 
(5) 
(6) 



with 

W e (k,r s 



= (-1)' 



4-7T 

21 + 1 



r 2 drlf f (r,r s )j,(/:r) (7) 



where the % has been inserted for later convenience. By 
making the filter compensated, we reduce the sensitiv- 
ity to the poorly constrained power at large scales and 
the dependence on the uncertain mean density in the 
sample. The corre lation functio n is defined such that 
t(r,fi) + 1 oc n' 2 (|PeeblesllT980l ). Integrating the left- 
hand side of this equation against We(r,r s )Le(fi)d 3 r re- 
sults in uJi{r s ) scaling directly with rT 2 following equa- 
tion Q. The constant term integrates to as the fil- 
ter is compensated. Hence, any uncertainty in n en- 
ters as a pure multiplicative offset in uje{r s ), which is 
less likely to overwhelm the acoustic signature at large 
scales. This in essence, eliminates sensitivity to the in- 
tegral constraint, which is a small effect to begin with. 
This feature of the filter also makes the statistic mea- 
sured in different subvolumes of a survey more indepen- 
dent. We expect that this will make internal error esti- 
mates f rom methods such as boots trap or jackknife more 

robust (IPadmanabhan et al.ll2009ri . 

Following IPadmanabhan etaU (12001 we consider a 
low order, smooth compensated filter. For simplicity we 
assume We is independent of £, though we could of course 
choose different weights for each multipolfl In terms of 
x = (r/r s ) 3 , the filter 



W(x) = (2a;) 2 (1 - xf 



1 



f 



(8) 



satisfies W(0) = W'(0) = W(l) = W'(l) = and 
J dx W(x) — 0. The suggested form in configuration 
space (top panel of Figure [1]) has a broad hump peak- 
ing at r ~ 0.65 r s that matches the width of the acous- 
tic peak in £(r) and a sharp (negative) spike at 0.9 r s 
of width 0(10%). This filter will smear a feature, such 
as the acoustic peak in £(r), by very little which means 
that the acoustic information will be localized in uje(r s ), 
however, not as localized as in £(r). Obviously, given 
sufficient signal-to-noise, measuring uie(r s ) for many r s 
values would allow resolution in £(r) even below the in- 
trinsic width of W(r/r s ) (see 



6 For example, we could make the fc-weight for £ = and 
1 = 1 equal. Doing so facilitates the computation of Q(r) = 

3/r *%J?}X r ,_t oM = ^0§$TT5)P flHamiltonl HOm when 
transformed to r a space, since it involves both the monopole and 
the quadrupolc. This ratio is useful for estimating the anisotropic 
parameter /3, however, it has limited applicability outside linear 
theory. 
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Fig. 1. — (top) The filter We(r/r B ) in configuration space. There 
is a broad hump matching the width of the acoustic peak in £(r) 
which peaks at r ~ 0.65 r s and a sharp (negative) spike at 0.9 r s 
with a width on the order of 10%. A filter of this shape will smear 
a feature, such as the acoustic peak in §(r), by only a small amount 
which means that the acoustic information will be well localized in 
LOeirs). Its compensated nature implies that LUi(r s ) is not sensitive 
to the integral constraint, (bottom) The filter Wi(kr s ) for £ = 0. 
The insensitivity of this filter to large scales is reflected in the 
fact that it is singly compensated and Wp ~ k 2 as k — > 0. At 

small scales We(kr s ) — > cos(fcr s )/(fcr s ) 4 . This is a much more 
rapid drop-off than observed in the kernel for ?(r), which scales 
as (fcr) -1 . These properties of the filter imply that o^(r s ) only 
probes a narrow range of scales in Fourier space, and that it is 
insensitive to large scale fluctuations or poorly constrained small- 
scale structure. 

With this choice of We the window function We can 
be computed analytically (see Appendix for numerical 
details) or numerically via fast Hankel transforms. We 

show We(kr s ) for I = 0, in the bottom panel of Figure 

[TJ Since the filter is singly compensated, We ~ k 2 as 
k —> 0, reflecting insensitivity to large scales. At small 
scales We(kr s ) — > cos(kr s )/(kr s ) 4: , a much more rapid 
convergence than evinced by the kernel for £(r), which 
scales as (kr)^ 1 . Thus uie probes a narrow range of scales 
in Fourier space and is insensitive to fluctuations on large 
scales or poorly measured or modeled small-scale struc- 
ture. One can choose the range of k to be sampled by 
appropriate choice of r s : more information from high k 
modes can be included by using smaller r s . 

As an example, the linear theory monopole statistic 
uo(f s ) is plotted in Figure [21 Plotting r 2 wo versus r s 
gives a convenient vertical range. The acoustic infor- 
mation is mostly localized into a single dip around the 
acoustic scale (see ^4.2.ip . 

2.2. Computation 
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Fig. 2. — Linear theory monopole statistic uo(r s )- The acoustic 
information can be seen around the acoustic scale, mainly localized 
into a single dip feature. 



It is possible to adapt ui£(r s ) into a sum over unbinned 
pair counts for any sa mple of mass tracers follow ing the 
methods described in iPadmanabha n et al.l (|2007l ); there 
is no need to first compute £ (r, ll) via binning of data. 
Pair counting allows us to record each galaxy's angle to 
the LOS (ll) accurately. Hence, like uJi(r s ) is sen- 
sitive to any anisotropics of the tracer distribution in 
clusters (discussed in SJTJ) . 

The redshift-space correlation function can be esti- 
mated as 

DD(r,n) 

&>( r )A») = "HB7 ( 9 ) 



RR(r,Li) 



where DD(r, ll) is the number of data tracer pairs sep- 
arated by r and have LOS angle corresonding to ll. 
i?i?(r, ll) is the analogue for randomly distributed points, 
normalized to the data counts by a factor of N^/N^. 
Here, Njj and Nr are the total number of data and ran- 
dom points respectively. When analyzing observational 
data, the number of random points needs to be much 
larger than the number of data points to keep the shot- 
noise in RR smaller than that in DD, especially at small 
r. For simulation data, however, it is not necessary to 
use a very large number of random points to compute 
RR smoothly at small scales (elicited below). 

Equation ([9]) implies that Equation ((4} can be rewrit- 
ten as 



Mrs) =i l / d a rWt(r)Lt(n) 



DDjr, li) 
RR(r,fi) 



(10) 



The —1 integrates to due to the compensated nature 
of the hlter. 

The RR piece is purely geometrical and is dependent 
only on the survey geometry (encoded in $(r, ll)) and the 
number of random points. Hence we can write RR as 



RR(r,Li) = 2Tm D N D r 2 $(r, Li)drd,Li 



(11) 



where no is the number density of data points which 
is easily calculable for surveys with well defined bound- 
aries. The above equation defines $(r, ll) to be any mis- 
match between infinite sized surveys/simulations and fi- 
nite sized ones due to the presense of boundaries. The 
noNjj factor is due to the normalization of the RR 
counts as mentioned above. For observations, <£>(r, ll) 
can be computed via binning methods and then fit using 



a smooth function $(r, li). Note that the binning men- 
tioned here is only required in the computation of the 
RR counts; there is no need to bin the data. In the case 
of simulations in a periodic box, $(r, ll) is constant as 
the volume is effectively infinite. 

With these points in mind, we can now pick arbitrarily 
small bins when computing the DD counts since RR has 
been approximated by a smooth function and hence does 
not suffer from shot-noise induced through pair-counting. 
As is such, we may employ a binning scheme in which 
there is either zero or one DD pair per bin. This step 
reduces the integral in Equation (|TU|) to a sum over DD 
pairs as in 



ut{r s ) 



iEDD 



Wj{ri)L(jM) 
n D N D V${ri,Lii) 



(12) 



Since the estimator can be written as a summation, there 
is no longer a need to bin data at all. 

2.3. Covariance Matrix 

Since uj£ does not require the binning of data, we can in 
principle estimate it at as many r s values as we wish with- 
out affecting the signal in the adjacent values: there is no 
bin which is made smaller. However adjacent points be- 
come increasingly correlated as the r s spacing decreases, 
compromising the usefulness of very fine sampling. 

In the Gaussian limit, the covariance matrix is 



Cov [u t (r s ),u t (r' a )] 
k 2 dk 



2(21 + 1) (21' + 1) 
V 



2tt 2 



Wi{kr s )Wv{kr l a )l w {k) 



(13) 



with 



Iw = 



dfj,Lt(fi)Le'(n) 



Y,Pdk)L L {n) 



(14) 

where -Pl(&)Ll(/z) is the legendre decomposition of 
the full 3D power spectrum P(k) and H is shot-noise. 
Assuming Poisson shot-noise, H = where n is the 

number density of the mass tracer. 

2.4. Summary of Key Features 

We conclude this section with a summary of the key 
features and advantages of uJi{r s ) over £(r) and P(k) 
estimators. 

1. u)g(r s ) has a compensated filter that reduces sensi- 
tivity to poorly constrained large scale power and 
hence the integral constraint. £(r), on the other 
hand, experiences these problems. The compen- 
sated filter also makes uji(t s ) measured in different 
subvolumes of the survey more independent which 
is important for attaining robust error estimates 
from methods such as bootstrap and jackknife. 

2. The filter is approximately compact in both con- 
figuration and Fourier space. The smoothness in 
configuration space leads to the steep drop-off at 
high k in Fourier space. This effectively minimizes 
the impact of large k or small scale power which is 
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not well constrained in large cosmological surveys. 
The filter is localized in configuration space which 
means that, unlike in P(k), the acoustic informa- 
tion is localized in ug(r s ). However, it is not as 
localized as in £(r) (see §4.2.1|) . 

3. Like u>e(r s ) can be easily adapted into a pair 
count statistic, so the relative positions of trac- 
ers and each pair's angle from the LOS is accu- 
rately recorded. Hence, it estimates the underlying 
galaxy distribution without the need to deconvolve 
a window function (as in the case of P(k)) and 
is sensitive to any anisotropics in the clustering of 
tracers. The anisotropics can be used to determine 
the underlying cosmology. P(k) estimators do not 
typically allow this type of analysis. 

4. There is no need to bin data when computing 
tog(r s ), unlike when estimating £(r) and P(k). 

While the setbacks of traditional P{k) and £(r) esti- 
mators are minor, it is beneficial to have a statistic like 
Wi{r a ) which combines many of the advantages of both. 

3. SIMULATIONS AND ANALYSIS METHODS 

A major goal in developing this new statistic is to use 
it for better measuring and calibrating the acoustic scale. 
The monopole statistic uJo(r s ) is especially useful in this 
regard as it is a direct map from the traditional 2-point 
correlation function and power spectrum. Therefore, any 
results obtained from BAO analysis via these three statis- 
tics is readily comparabl e. As with the quadrupole of th e 
power spectrum P 2 (k) ((Padmanabhan fc White] 12008(1 . 
anisotropic BAO analysis can be performed using the 
quadrupole u>2(r s ). However, as the main goal of this pa- 
per is to give a broad introduction to the u)e(r s ) statistic 
and a simple demonstration of its application, we defer 
detailed discussion of oj^i^s) and anistropic BAO to a 
future paper. 

In this section, we implement uJo(r s ) through the use 
of pure dark matter N-body simulations. It is possi- 
ble to model a variety of different galaxy populations 
and biases through application of appropriate HODs to 
the halos found in the simulations. This is important in 
demonstrating the robustness of the ujt{r s ) statistic over 
a diverse set of galaxy populations. We compute uio(r s ) 
for dark matter and the mock galaxy populations created 
via the HODs, and demonstrate how it can be used to 
measure the shift in the acoustic peak. This quantity is 
important in constraining the precise size of the acoustic 
scale, which may be slightly different from that predicted 
by linear theory due to nonlinear structure growth. In 
order to quote the shift with accurate errors, we use a 
resampling technique described in §3.41 which gives us a 
large number of shifts from which to calculate a mean 
and a standard error of that mean. 

3.1. Simulations 

Our simulations were performed using a code devel- 
oped by Metchnik and Pinto w hich employs a new 
method ((Metchnik fc Pinto! f20Toh to compute N-body 
forces under periodic boundary conditions. Rather than 
resorting to the Fourier methods of PM schemes or using 
Ewald sums to represent periodicity, this new method 



represents periodicity directly, expressing the force on a 
particle as due to the rest of the simulation volume and 
an infinite sum over its periodic images. 

This is made more efficient by partitioning the com- 
putational domain into a three-dimensional grid. The 
acceleration on particles within a grid cell is divided into 
two parts: a near field and a far field. The near field is 
that due to the other particles in the cell and to particles 
in the adjacent 26 cells. In these calculations, the near- 
field acceleration was computed using the direct, 0(N 2 ) 
method, with Plummcr softening. 

The far field acceleration on particles in the cell due 
to each more distant cell on the grid is represented as a 
Taylor series expansion based on the multipole moments 
in the distant cell. The contribution from all periodic 
images of the distant cell is included by recognizing that 
the multipole moments in a cell are identical to those in 
all of its images. Thus, the sum over images depends 
only on the (fixed) geometry of the grid and need be 
performed only once. The (relatively) small set of val- 
ues which results provides a simple and rapidly-evaluated 
relation between the multipole moments in distant grid 
cells and the Taylor coefficients of the expansion for the 
acceleration in a given cell, all under periodic boundary 
conditions. The acceleration due to all of the periodic 
images of the cell and its 26 neighbors is included in a 
similar manner. 

The calculations described here used order-16 expan- 
sions, providing an overall force accuracy per particle 
which agrees with Ewald summation to better than six 
decimal digits. Overall, the method is significantly faster 
and requires significantly less memory than other meth- 
ods for computing accelerations in large, periodic N-body 
problems, while providing accelerations to machine pre- 
cision for all particles. Time integration was performed 
using a standard kick-drift-kick algorithm, using inde- 
pendent time steps for each particle. 

We derive our results in redshift-spacc at z = 1 from a 
set of 44 simulations with 1024 3 dark matter particles in 
each. The simulations were performed in l/i -1 Gpc peri- 
odic boxes with the WMAP5+SN+BAO best-fit cosmo- 
logical parameters: O m = 0.279, il\ = 0. 721, h = 0.701, 
Q b = 0.0462, n s = 0.96 and cr 8 = 0.817 (|Komatsu et al.l 
12009(1 which implies a particle mass of 7.2 x 10 w h- 1 M Q . 
The initial conditions are generated via t he second-or der 
Lagrangian perturbation theory code of iSirkol (|2005H at 
z = 50 with no extra power for the box scale. 

3.2. Halo Occupation Distributions 

We use a simple friends-of- friends ((Davis et alj [l985) 
method with a linking length equal to 0.16 of the in- 
tcrparticle spacing to identify the collapsed dark matter 
halos in our simulations. We then populate these ha- 
los with galaxies by applying simple HODs based on the 
form 

(N g (M)) = [1 + (M/M sat )T]exp(-M cen /M) (15) 

where M is the halo mass, M ccn is the minimum 
mass for a halo to contain a central galaxy, M sat 
is the minimum mass for a halo to contain at least 
one satellite, and 7 is an exponential parameter that 
we set to 1 dGuzik fc Seliakl 12001 iBerlind et aL| [200l 
iKravtsov et al.l 120041: iZheng et al.l I2005H . We assign a 
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central galaxy to a halo if M > Af ccn , this is a good ap- 
proximation to (N ccn (M)) = exp(— M cen /M). The cen- 
tral galaxy is taken to be at the halo's center of mass 
and assigned the center of mass velocity. If a halo is 
assigned a central galaxy, then the number of satellite 
galaxies is determined by generation of a random in- 
teger based on a Poisson distribution with mean equal 
to (iV sat (M)) = (M/M sat ) 7 . We then randomly pick a 
corresponding number of halo particles and assign their 
positions and velocities to the satellites. 

In order to compare the peak shifts derived using 
u)o(r s ) and P(k) in a range of models, we apply three 
different HODs to our simulations. The properties of 
these HODs are described in Table [TJ and obtained by 
adjusting the values of M cen and M sat . We also list 
the properties of a dark matter (DM) only case that is 
merely a 0.4% subsample of the particles in each sim- 
ulation. We compute u>o(r s ) in redshift-space via the 
pair counting method detailed in £12.21 within the range 
bh,- 1 Mpc < r a < 200/i _1 Mpc using 5/i _1 Mpc spacing. 
We also compute as, the RMS mass fluctuation within a 
8/i _1 Mpc radius, using a similar pair counting method 
derived from t he configurat i on spa ce equation for a gen- 
eral radius R (jZehavi et all (|2005t )) 
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3.3. Fitting wo(r s ) to Measure the Peak Shift 
We fit the redshift-space LJo(r s ) using the form 



P sim (k) = B{k)P m (k/a) + A(k) 



where 



B(k) 



(6| + b\k + b 2 k 2 ) 

1 ~t~ 7 scaled 



(17) 



(18) 



and A(k) transforms into A(r s ) = a\rj 9 in r s space. 
Such a form for A(k) is motivated by the fact that 
we want to marginalize over the shape of the correla- 
tion function at small scales (i.e., the contribution of 
the 1-halo term). Expanding Equation ((HJ, we see that 
W(x) = (-4x 5 + 2.t 4 - Sx 3 + 2x 2 )^. This implies that 

the highest order term in r s is r~ 9 . At small scales, this 
is the term that will dominate in the transformation from 
£(r) to uJo(r s ) as defined by Equation (|4]). 

P m {k) is the template power spectrum we use for our 
fitting. To account for the degradation of the acous- 
tic peak through nonlinear evolution and redshift distor- 
tions, the template model P m (k) is obtained from the 
linear power spectrum Pn n (k) at z = 1 by the modifica- 
tion 

Pn{k) = [Fi in (fc)-P sm0 oth(fc)]exp(-fc 2 I]2 1 /2)+P smooth (fc) 

(19) 

where P S mooth is the dewiggled power spectrum de- 
scribed in lEisenstein &: Hul (|1998D and S n i is a nonlinear 
parameter used to degrade the peak (lEisenstein et all 
2001 iTegmark et all [2001 IC rocce k, Scoccimarro! 12006; 
Eisenstein et all 120071 ; ICrocce fc Scoccimarro! 120081 ; 
Matsubaral!2'008l) . To allow maximum flexibility in our 



marginalization, we marginalize over E„i and the other 
nonlinear nuisance parameter in equation (|18[) . r sca i c . 

Our scale dilation parameter a represents the shift in 
the acoustic peak. Under this formalism, a > 1 indicates 



a shift towards smaller scales and a < 1 indicates a shift 
towards larger scales. Physically, a is the ratio between 
the linear theory acoustic scale (150 Mpc) to the mea- 
sured acoustic scale. Since all the terms in the fitting 
function are additive, the basis functions P m (fc/a)/(l + 
fscaicfc), kP m (k/a)/(l + r sca \ c k), k 2 P m (k/a)/(l+r sca \ c k), 
and r~ 9 for a fixed a, £ n i and r sca io can be easily mapped 
into r s space (if necessary) using equation ([6]) . A least- 
squares fit using the mapped basis functions is then per- 
formed against ujq{t s ) from the simulations to obtain val- 
ues for the linear nuisance parameters b\, 62, &3 and a\. 
As we are interested in the acoustic feature, we use a fit- 
ting range of 30 < r s < 200/i _1 Mpc. For an r s spacing 
of 5/i _1 Mpc, this implies 28 degrees of freedom in the 
fit, where the number of degrees of freedom is defined as 
the difference between the number of data points being 
fit and the number of parameters in the fitting form. 

We assume that the errors on uio(r s ) can be well ap- 
proximated by the covariance matrix C assuming Poisson 
shot-noise (see ^3.4|) with the addition of nonlinear shot- 
noise (see equation (|20|) and surrounding text). We also 

assume that the monopole (£ = 0) dominates P(k) so 
that all higher order contributions to the power spectrum 
are effectively zero. This amounts to computing C using 
P(k) = P m (k) as the input power spectrum, where we 
take a fix ed E n i = 7 . 0fe~ 1 M pc at z = 1 in redshift-space 
following ISeo et all (|2008|) . We normalize this power 
spectrum to the amplitude of the redshift-space power 
spectrum through multiplication by the bias squared de- 
fined initially as b 2 = (o-s.casc/cs.mattcr) 2 - The values 
for 

"8. case are given in Table [2] and (7s tm atter = 0.506 in 
real space at z — 1 in linear theory. We want to ensure 
that the input power spectrum to the covariance matrix 
calculation is as close as possible to the simulation data 
so that the covariance matrix is a reliable estimate of 
the errors. To do this, we marginalize over the average 
wo(r s ) of the 44 simulations for each case to obtain a 
value for the leading order term that scales P m (k): b 2 in 
the fitting form of equation (|T5|) . We then iterate this 
marginalization and scale b 2 by the values of b\ obtained 
until the output b\ from this iterative fitting is close to 
1. We expect that scaling P m {k) by this final value of b 2 
will approximate the simulation data well and hence be 
valid input to the covariance matrix calculation for the 
resampling techniques described in £13.41 

The shot-noise we enter into the calculation of C in- 
cludes a nonlinear component (quoted in Table [T]) in 
addition to the Poisson shot-noise n _1 as described in 
Equation (|T3"|) . This additional shot-noise is a result of 
nonlinear structure formation on small scales. We esti- 
mate this nonlinear shot-noise as 



4^r 2 [£(r)-6in(r)]dr (20) 



where £(r) is the correlation function averaged over the 
44 simulations for each HOD, £iin(0 is the linear corre- 
lation function at z = 1 and r 110 niin is the scale above 
which nonlinear effects become unimportant. We take 
rnoniin to be 10ft" 1 Mpc. The resulting K non ij n is a rough 
estimate of the excess small scale correlation due to non- 
linear evolution. Since it makes little difference whether 
all of this extra shot-noise comes in at zero separation in 
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Fig. 3. — uo( r s) averaged over all 44 simulations for HOD1 (black 
diamonds). Overplotted are a fit obtained through the form in 
Equation H17H (red crosses) and a Oth order fit in P m (k/a) (purple 
dots), both over a range of 30 < r s < 200h~ 1 Mpc. One can see 
that the Oth order fit already appears quite good with \ 2 = 2.20 
per dof. However, by introducing additional nuisance parameters, 
the quality of the fit over the specified range improves further to 
X 2 = 0.94 per dof. 

r or through the extended effects of the one-halo term, 
which is only important at small r, we assume the excess 
correlation to be a spike at r = for convenience. When 
transformed into k space, this gives a constant and be- 
comes extra white noise that we add on to every mode 
equally, in addition to the Poisson shot-noise. For the 
subsampled DM case, the linear and the measured corre- 
lation functions were sufficiently similar at 1-lOft -1 Mpc 
to warrant taking K non ii n = for this case. Alternatively, 
one can also account for nonlinear shot-noise by com- 
puting the covariance matrix using the nonlinear power 
spectrum, but this is more computationally challenging. 

The x 2 likelihood indicator corresponding to the best- 
fit linear nuisance parameters for fixed a, £ n i and r sca i c 
is then 

X 2 = (LUo-mfC-'iujo-rn) (21) 

where ujq is u)o(r s ) measured from the simulations, rn is 
the best- fit model and C _1 is the inverse of the covari- 
ance matrix. We compute the best-fit values of a, E n ] 
and r sca i e by minimizing \ 2 of the fits for the DM case 
and for each HOD using a generalized reduced gradient 
method from IDL. We quote the bias for each case as 
C8,casc/o8, matter multiplied by the additional scaling fac- 
tors of b\ described above. 

We plot wo(fs) averaged over all 44 simulations for 
HOD1 in Figure [3J Overplotted are the marginalization 
obtained through the form in Equation (|17l) and a Oth 
order fit obtained by just a rescaling of P m (k/a), i.e. 
b 2 P m (k/a) where b is the only fit parameter. Although 
the 0th order fit already appears quite good, one can 
see that introducing additional nuisance parameters im- 
proves the quality of the fit even more over the fitting 
range. The \ 2 P er degree of freedom (dof) improves from 
2.20 to 0.94. 

Using this fitting technique, we derive values of bias, 
a and associated errors for each of our three HODs and 
our DM case via the resampling methods described in 
the following section. 

3.4. Resampling Methods 



We use two different methods to measure the mean 
peak shift a and the scatter in the mean a a for each case 
in Table [TJ The first is a modified jackknife technique 
in which we randomly select M out of N simulations at 
a time without replacement, average their ujo(r s ) and fit 
this average. We repeat this 1000 times and extract an 
average a and a scatter in a. This scatter needs to be 
rescaled by an additional factor of / = \J~hI / \J N — M in 
order to reflect the scatter asscociated with the mean of 
a for N simulations. For our simulations we have N = 44 
and take M — 22. With this choice of M, f — 1 and so 
the scatter in a reflects the error in the mean of a. This 
method is useful in that it provides us with a large set of 
a's from which we can accurately derive a mean a and 
<r a . 

The fit results for the subsampled DM CcLSC clS well 
as for each HOD model are quoted in Table [2] The 
average values of £ 1U and r sca i e are also included for 
completeness, however the focus of this paper is on 
a. The values of a we obtain are la consistent with 
those derived using the pertur bation theory results of 
iPadmanabhan fc White] (|2009f ). The somewhat low 
value of x 2 P er dof for HOD3 suggests that we are over- 
estimating the amount of nonlinear shot-noise. We are 
also approaching the shot-noise limited regime for HOD3 
as evidenced by the fact nPo.2 ~ 1.6. We note here that 
we used % 2 only to find the best-fit a for each HOD, not 
to generate the errors. Hence, the fact that our reduced 
\ 2 values are slightly deviant from unity does not hinder 
the error estimation. 

The second method we use is jackknife resampling. 
The results obtained using this method are in good agree- 
ment with the first method. This indicates that the error 
estimates obtained using our first method are robust in 
comparison to more traditional methods. 

By using these resampling techniques, any non- 
Gaussian effects not accounted for by assuming a Gaus- 
sian covariance matrix while fitting (as in £|3.3p will be 
reflected in a a . 

4. COMPARISON TO THE POWER SPECTRUM 
4.1. Comparison from Simulations 

An important step in implementing this new statis- 
tic is to show that it produces consistent results when 
compared to established methods and can therefore be 
an effective calibrator of the acoustic scale. We do this 
by comparing the peak shifts measured from the same 
set of simulations via the new uio(r s ) statistic and the 
traditional P{k) method (see S&M for details). 

As we use the same 1000 random sets of M simulations 
as S&M, there should be a 1:1 correspondance between 
the a's derived from uJo(r s ) and P(k) for the DM case 
and for each HOD. It should be noted that S&M use dif- 
ferent P(k) fitting forms from the one detailed in i)3.3l 
They employ two fitting forms, both of which can also be 
described by equation (|17[) . The first form has B(k) as a 
2 nd order polynomial and A(k) as a 7 th order polynomial. 
The second form uses Pade approximants for B(k), i.e., 
B{k) = 6 (l+cifc+C3fc 2 +C5fc 3 )/(l+c 2 fc+C4fc 2 ) and a 2 nd 
order polynomial for A{k). We have chosen a different 
form in this work to induce better convergence of the in- 
tegral from Fourier space to r s space while transforming 
the basis functions. We compare the a's measured from 
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Fig. 4. — a from Loo(r B ) versus a from P(k) for HOD2 (top) and 
DM (bottom). The data points are from a resampling technique 
in which we randomly pick M = 22 simulations out of N = 44 
total and fit the averaged uJo(r s ) from these M simulations. We 
repeat this 1000 times and hence obtain 1000 values of a. The 
scatter on a needs to be rescaled by yM / y/N — M to reflect the 
true scatter on the mean. For our choice of M, this scaling factor 
is equal to 1. Hence the scatter in the plot truly reflects the scatter 
on the mean of a. The red cross marks the mean a values with 
their associated errors. The central grey line has unity slope and 
passes through the mean. The two outer grey lines delineate the lcr 
boundaries associated with Aa^p- As the data points lie largely in 
between the la lines with a slope similar to unity for both HOD2 
and the DM case (see Table [3}, we conclude that the two a sets 
are consistent with each other. The same correlation is observed 
for HOD1. HOD3 shows 1.6rr agreement between u>o(r s ) and P(k). 
This slightly larger discrepancy may be due to the fact that shot- 
noise is becoming significant in this low number density case. Also, 
shot-noise may affect P(k) and LJo(r 3 ) differently or the nuisance 
parameters may not be fully handling the scale-dependence of a 
high-bias HOD such as HOD3. The large scatter in the DM case 
is likely due to the subsampling of matter in the computation of 
^o(^s) but not in P(k). The overall agreement between the uio(r s ) 
and P(k) results imply that distance measures will be consistent 
between the two. 



uig(r s ) against th ose measured from P( k) by ISeo et ahl 
(|2010f) (DM) and iMehta et all d in preuD (HODs) using 
the first form. It should also be not ed here that the P(k) 
results obtained for the DM case bv lSeo et al.l (l2010l) uti- 
lize the full DM sample whereas we have subsamplcd to 



reduce computation time in this work. The methodology 
used to derive th e P(k) results are described in detail in 
ISeo et all pOlOh . 

Figure 0] shows a from uio(r s ) versus a from P(k) for 
the 1000 fit iterations performed on HOD2 (top) and 
DM (bottom). The red cross indicates the mean a val- 
ues with their associated errors. The central grey line 
is a line with slope unity that passes through the mean. 
The two outer grey lines indicate the la boundaries as- 
sociated with Aa u p = a u — ap. The mean difference 
between a u and ap along with the standard deviation 
is quoted in Table [3] under (Aa^p) for the DM case as 
well as each HOD case. The plots indicate that the cor- 
relation between a u and ap is la consistent with a line 
of slope unity that has a y-intercept of for both HOD2 
and the DM case. This implies that the two a sets are 
consistent with each other. The same holds for HOD1 as 
evidenced through the table, but HOD3 is slightly more 
deviant with 1.6a agreement between wo( r s) and P(k). 
The larger discrepancy between the HOD3 results may 
be due to the fact that shot-noise is becoming significant 
in this case (as shown in fc|3.4[) . It could also be that shot- 
noise affects P(k) and uJa(r s ) differently or the nuisance 
parameters are not fully handling the scale-dependence 
of a high-bias HOD such as HOD3. The large scatter in 
the DM case is likely due to the fact that we have sub- 
sampled the matter in our computation of wo(r s ) but not 
in P(k). The en's from most of the cases are la consistent 
between wo(f s ) and P(k), indicating that distance mea- 
sures will be consistent between the two statistics. This 
also indicates that any systematics introduced by using 
the different fitting forms for P(k) and uio(r s ) are minor. 
Hence we conclude that wo( r s) is a well-tuned statistic 
for analysis of BAOs. 

4.2. Theory Constraints on a a 

As we wish to promote uo(r s ) as an alternative method 
for analyzing the BAO, it is necessary to show how much 
acoustic information can be extracted from wn(r s ) rela- 
tive to P(k) and £(r) for surveys of the same size. If 
our fitting ranges were infinite, then by the definitions in 
iJ2J all three estimators should yield the same amount of 
BAO information. However, in reality, fitting ranges are 
finite. 

We investigate the effects of this by shifting the acous- 
tic feature in the linear theory uj (r s ), P(k) and £(r) by 
a given a and then running our fit algorithms to see how 
well we can recover this input a. For wo(r s ) we fit be- 
tween r s = 30-210ft -1 Mpc in r s spacings of 2.hh~ 1 Mpc; 
for P(k) we fit between k = 0.0-1. 2h Mpc" 1 in log(fc) 
spacings of ~ 0.002; and for £(r) we fit between r = 20- 
200ft, -1 Mpc in r spacings of 1ft -1 Mpc. We use the same 
fitting technique as described in i j3.3l but with different 
forms for B(k) and A(k) that are then transformed to 
r and r s space to fit £(r) and wo(r s ) respectively. This 
means that the fitting forms for P(k), £(r) and u>o(r s ) all 
derive from the same B(k) and A(k) functions. By en- 
forcing this consistency between fitting forms, any poten- 
tial systematics that may arise due to the use of different 
fitting forms for each statistic can be avoided. 

We pick B(k) = b, where b is the large-scale bias, 
and A{k) to be the cold dark matter-only power spec- 
trum multiplied by a set of cubic spline functions spec- 
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Fig. 5. — x 2 versus a for fits in LJo(r s ), P(k) and £(r). We 
shift the acoustic feature from linear theory by a given a (here 
Cgiven = 1-0) and run our fit algorithms to see how well we can 
recover this input. For u;o(r s ) we fit between r s = 30-210ft _1 Mpc, 
for P(k) we fit between k = 0.0-1. 2ft Mpc -1 and for £(r) we fit 
between r = 20-200ft -1 Mpc. The parabolic shape of the curves 
is due to the fact that uio(r s ), P(k) and £(r) are derived from a 
Gaussian random field in linear theory. The width of the parabola 
at x 2 = 1 is then the theoretical a a of the fit. The uJo(r a ) and 
P(k) curves overlap nicely, implying that the o a ratio between 
bjQ(r g ) and P(k) is ~ 1. This indicates that for the given finite 
fitting ranges, u)o(r s ) and P(k) contain equal amounts of acoustic 
information (trivially true for infinite fitting ranges by definition of 
o;o(r s ), however, real surveys are finite in size). Hence, the same 
amount of acoustic information can be obtained through either 
<^o(fs) or P(k) analysis of equal volume surveys. 

ified at k = 0. 01-1. 09/i Mpc -1 . The spline points are 
picked so that seven of them are logarithmically spaced 
in the range k = 0.01-0.25 and seven of them are linearly 
spaced in the range k = 0.25-1.09 giving a total of 14 
spline points. The basis functions specified by the small 
k spline points are necessary to allow flexibility in the 
marginalization of u>o(r s ) at large scales. The derivative 
is also specified at the first spline point to derive an ad- 
ditional spline function. The spline functions are taken 
to be natural (i.e., second derivative equal to 0) at the 
first and last spline points, beyond which linear extrap- 
olations are implemented. This choice of A{k) ensures 
convergence when transformed to r s space and makes 
the fits in P(k), £(r) and u>o(r s ) readily comparable. We 
assume a survey volume of l/i -1 Gpc with H = 1000 
(i.e., one million particles). As a cross check, we con- 
firmed that this new fitting form docs in fact give similar 
results to the form used in .31 

Figure [5] plots the x 2 versus a for wo(r s ), P(k) and 
£(r). Here the input a is equal to 1. If w (r s ), P(k) 
and £(r) are derived from a Gaussian random field as 
they are in linear theory, we would expect that x 2 ver- 
sus a be parabolic as shown in the figure. The width 
of the parabola at \ 2 — 1 is then the theoretical a a 
from the fit. The overlap between the coo(r s ) and P(k) 
curves indicates that the ratio of a a for ojq{t s ) to P{k) is 
~ 1. This means that for the given fitting ranges u>o(r s ) 
and P(k) contain equal amounts of acoustic information. 
Since volume is proportional to \ 2 which is proprotional 
to er~ 2 , an important implication is that we are able to 
obtain the same amount of acoustic information using ci- 
ther uio(r s ) or P(k) analysis of equal volume surveys. We 
emphasize here that the results presented in this section 
assume idealized linear theory forms for P(k), £(r) and 



u>o(r s ). In practice, numerous physical and observational 
effects distort the measured statistics from these ideals. 
However, we expect that the features described in ^2] will 
reduce the impact of troublesome observational effects in 
any BAO analysis using wo(r s ). 

4.2.1. Locating the Acoustic Information 

It is useful to track down where the acoustic informa- 
tion lies and how it changes with a in u>o(r s ), and 
P{k). This is reflected in the derivatives duo/da, d£/da 
and dP/da after marginalizing out the broadband shape. 
To do this, we calculate the residuals from the P(k) fits 
described in JO] for a = 0.996 and a = 1.004. These 
residuals should be representative of the acoustic signa- 
ture after the broadband shape has been marginalized 
out. We then take dP/da as the difference between these 
residuals divided by 1.004 - 0.996 = 0.008. The trans- 
formations of dP/da into r s and r space then give us 
dujo/da and d^/da respectively. We have plotted dP/da 
in the top panel of Figure El d^/da in the middle panel 
and dtdo /da in the bottom panel. If one plots the ratio of 
dP/da to P(k)/k, one is left with the approximate shape 
of the signal-to-noise ratioQ (SNR). This is plotted in the 
top right-hand corner of the top panel in Figure [BJ One 
can see that the SNR is small at k < 0.05/i Mpc -1 , indi- 
cating that the small k ringing in dP j dot is merely noise 
from the spline basis functions attempting to match the 
shape of P(k) at these scales, and is not indicative of 
the shift in acoustic information with a. The shifting 
of the acoustic information with a is only truly evident 
at k > 0.05/i Mpc -1 where the SNR is larger. It is ev- 
ident from these plots that the acoustic information is 
not as localized in uio(r s ) as in £(r), but it is still rea- 
sonably well localized. The bottom panel of Figure [S] 
indicates that the optimal fitting range that will include 
all of the acoustic information encoded in u>o(r s ) is some- 
where within the range r s = 30-300/i _1 Mpc. 

The top panel of Figure [7] shows how a a changes as the 
minimum r s of the fitting range is stepped up from 30- 
80/) . Mpc with the maximum r s of the fitting range 
fixed at 300/i" 1 Mpc. The bottom panel of Figure 
shows how a a changes as the maximum r s of the fit- 
ting range is stepped down from 300-120/j -1 Mpc with 
the minimum r s of the fitting range fixed at 30/i _1 Mpc. 
The regions of the plots where a a begins to modulate are 
blown up for clarity. One can see that when the mini- 
mum of the fitting range is larger than 40/i _1 Mpc, a a 
begins to deviate, indicative of missing some of the acous- 
tic information. This also happens when the maximum 
of the fitting range is smaller than 200/i _1 Mpc. Hence, 
the minimum fitting range that allows one to extract 
all of the acoustic information appears to be r s = 40- 
200/j _1 Mpc. Hence, to ensure we are encapsulating all 
of the acoustic information, we picked the fitting range 
to be r s = 30-210/i -1 Mpc in Figure 03 

5. CONCLUSIONS 

We have presented a new statistic uJe(r s ) for analyz- 
ing baryon acoustic oscillations. This new statistic is 

7 The noise term ap = P(k) / V dN where dN = k 2 dk is the 
number of modes out to k. For constant increments in fc, dN oc k 2 
and hence ap oc P(k)/k 
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Fig. 6. — (top) dP/da calculated from the residuals of the P(k) 
fits after marginalizing out the broadband shape. The change in 
P(k) with a captured by dP/da should correspond to how the 
acoustic information is shifted as a changes. The ratio of dP/da 
to P(k)/k approximates the shape of the signal-to-noise ratio and 
is shown in the top right-hand corner of the plot. The fact that 
the ratio is very small at k < 0.05h Mpc -1 indicates that all the 
ringing in dP/da at these scales is just noise from the spline ba- 
sis functions attempting to match the broadband shape of P(k) 
at these scales. Hence, this small k ringing does not actually re- 
flect the shifting of acoustic information as a changes. The ratio 
is larger for k > 0.05/i Mpc -1 indicating that the oscillations in 
dP/da at these k truly reflect the shifting of acoustic information 
with a. (middle) d^/da obtained by transforming dP/da. This 
shows where the acoustic information is located in configuration 
space and how it changes with a. (bottom) duio/da obtained by 
transforming dP/da to r s space. This shows where the acous- 
tic information is located in r s space and how it changes with a. 
Comparison with the middle panel indicates that the acoustic in- 
formation is not as localized in cuo(r a ) as it is in £(r), however, it 
is still reasonably well localized. All of the acoustic information is 
located within r s ~ 30-300/i _1 Mpc, indicating that the optimal 
fitting range for ujoi^s) is somewhere within these limits. 
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Fig. 7. — (top) Plots how cr a changes as the minimum r s of the 
fitting range is stepped up from 30-80/i — 1 Mpc with the maximum 
r s of the fitting range fixed at 300/i _1 Mpc. The deviation of a a at 
fs,min larger than 40/i _1 Mpc indicates that some of the acoustic 
information is being missed by these fitting ranges, (bottom) Plots 
how a a changes as the maximum r B of the fitting range is stepped 
down from 300-120/i — 1 Mpc with the minimum r s of the fitting 
range fixed at 30/i Mpc. The deviation of a a at r s , maa; smaller 
than 200h~ 1 Mpc indicates that some of the acoustic information 
is being missed by these fitting ranges. The above analysis implies 
that in order to extract all of the acoustic information, one needs 
to fit between r s = 40-200/i _1 Mpc at minimum. 



advantageous over the traditional methods used to es- 
timate £(Y) and P{k) as it does away with many of 
their setbacks. Estimators of £(r) are sensitive to poorly 
measured large scale power through effects such as the 
integral constraint, whereas the compensated nature of 
the filter W^(r) used to compute u>e(r s ) circumvents this 
problem. We expect that this feature of the filter will 
also make 0Je(r s ) measured in different subvolumes of a 
survey more independent which makes error estimation 
methods such as bootstrap and jackknife more robust. 
Estimators of P(k) give the true density field convolved 
with a window function making the measured statistic 
biased. Attempting to deconvolve the window function 
introduces artificial smoothing at small separations in k 
due to their near delta function shapes in real obser- 
vational surveys. u^(r s ), on the other hand, does not 
suffer from this as it is a pair count statistic. Pair count- 
ing estimators allow us to record the relative locations 
of tracer pairs as well as their angles from the LOS di- 
rection. This provides us with an accurate map of the 
underlying galaxy distribution. The fact that LOS an- 
gle is recorded also means that any anisotropic clustering 
should be apparent in u!((r s ). Hence we expect that it 
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can also be used to probe the underlying cosmology. In 
addition, there is less need to worry about binning re- 
lated issues when computing u>e(r s ) as we never need to 
bin the data. The smoothness of the filter in configura- 
tion space causes the rapid fall-off of the filter in Fourier 
space. This reduces the impact of large k modes or small 
scales which are not well constrained in large cosmology 
surveys. The localized nature of Wg(r) is conducive to 
minimal smearing of the acoustic information so that it 
is mostly concentrated in a single dip around the acous- 
tic scale. This translates to a cleaner representation of 
the acoustic information when plotted, as opposed to the 
oscillatory features of P(k). 

We also showed that with the present form for Wg(r) 
and a finite fitting range encompassing the acoustic scale, 
it is possible to extract the same amount of acoustic in- 
formation using either uo(r s ) or P(k) from equal vol- 
ume surveys. It is important to note that these results 
were obtained through analysis of idealized linear theory 
forms for P(k), £(r) and tJo(r s ). In practice, the mea- 
sured forms of these statistics are distorted by various 
physical and observational effects. However, we expect 
that the features described in $2] will reduce the impact 
of troublesome observational effects in any BAO analysis 
using ujQ(r s ). We also demonstrated where the acoustic 
information is located in wo(r s ), £(r) and P(k) and how it 



changes with a. From this analysis, the minimum fitting 
range required to extract all of the acoustic information 
from u>o(r s ) appears to be r s = 40-200/i _1 Mpc. 

We compared the acoustic peak shifts derived using 
u)o(r s ) to those derived using P{k) for a pure DM case as 
well as for three halo based galaxy models. The results 
for the DM and the higher number density cases are all 
much better than la consistent with each other. The low 
number density case is slightly deviant with 1.6a agree- 
ment between ujo(r s ) and P{k). This may be a result of 
approaching the shot-noise limited regime or our lack of 
understanding of shot-noise in general and how it may 
affect P(k) and uio(r s ) differently. It may also be caused 
by our fitting form not handling the scale-dependence of 
high-bias models in full. The general consistency betwen 
u>a(r s ) and P(k) is encouraging and implies that distance 
measures will be consistent between the two methods. 
From this and the features listed above, we conclude that 
ujo(r s ) is a well-tuned new statistic for BAO analysis. 
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TABLE 1 

HOD PROPERTIES 



Model 


Total # 
of Galaxies^ 


Satellite 
Fraction (%) 


(ft -1 M ) 


M sat 

(A _1 M ) 


"1 

(h 3 Mpc" 3 ) 


nPo.2 


^nonlin 


DM 


4 x 10 6 








0.004 


4.60 


0.0 


HOD1 


2 x 10 6 


5 


1.4 x 10 12 


9.2 x 10 13 


0.002 


5.78 


450.0 


HOD2 


1 x 10 6 


5 


2.6 x 10 12 


1.5 x 10 14 


0.001 


3.59 


700.0 


HOD3 


3 x 10 5 


5 


6.4 x 10 12 


3.1 x 10 14 


0.0003 


1.59 


1550.0 



Note. — HODs are referred to by the designations under the "Model" heading throughout the 
paper. 



a Number of DM particles in the DM only case. 
b The nominal Poisson shot-noise is n _1 . 



TABLE 2 
Fit results for each HOD model 



Model 


a - 1 






i scale 


x 2 


bias a as 




(%) 


(%) 






(per d.o.f) 


(bf 


DM 


0.0457 


0.2333 


6.66 


19.99 


0.92 


1.25 0.63 


HOD1 


0.1065 


0.2243 


5.61 


19.78 


0.94 


2.04 1.11 


HOD2 


0.1634 


0.2449 


5.85 


19.98 


0.86 


2.28 1.25 


HOD3 


0.4897 


0.3326 


6.27 


20.04 


0.72 


2.77 1.55 


Note. 


— Fitting ranj 


;e: 30 < r s 


< 200k- 1 


Mpc. a a 


is the error on the 


mean a of the 



44 simulations. 

a Bias is not equal to 1 for the DM only case because we are working in redshift space. 



New Statistic for Analyzing BAO 



13 



TABLE 3 
Difference in mean a between 
uj (r s ) AND P(k) 



Model (Aa„ P ) 
(%) 



DM 


0.0516 


± 


0.1205 


HOD1 


0.0076 


± 


0.0672 


HOD2 


0.0205 


± 


0.0600 


HOD3 


0.1035 


± 


0.0665 



APPENDIX 
EVALUATING W e (K) 

The expressions for Wg(k), in terms of polynomials of k times trigonometric functions, involve a lot of cancellation. 
This makes them unstable to direct evaluation. However if we define 



then 



while 



K n (k) 



W (k) = y 



2 + n 



(kr s 



a+n 



A:r s 



' sin x dx 



- (K 7 - K 16 ) - (K w - K 13 ) 



24n 

W 2 (k) = —^2 [3K 5 - 16K 8 + 25K n - 12K U ] 



It is straightforward to evaluate K n (x), the limits are 

n + 2 



KJx) = 1 



2 , n+2 4 
■x H — —, tX ■ 



3!(n + 4) 5!(n + 6) 



as x — > and 



: . . cosx , „. since 

K n {x) = -(n + 2)— 5- + n(n + 2)^- + 



as x — y oo. The K n also satisfy a simple recurrence relation 



K n (x) 



n + 2 



[nsina; — xcosx — (n — l)xK n -2\ ■ 



(Al) 
(A2) 
(A3) 

(A4) 
(A5) 

(A6) 



Use of this recurrence relation for high k and the power-series expansion for low k results in stable evaluation of the 
W e . 



